home *** CD-ROM | disk | FTP | other *** search
/ Ian & Stuart's Australian Mac 1993 September / September 93.iso / Archives / Applications / Calculators / NumberCrunch 1.41 / Number Crunch II Demo / Library Routines / Text Examples / Interpolation, Min, Zeros < prev    next >
Encoding:
Text File  |  1992-09-03  |  7.6 KB  |  234 lines  |  [TEXT/NCII]
  1. ################
  2. # Interpolation, Minima, and Zeros.
  3. # 
  4. #
  5. #   function Interpolate(X,Y, Xinterp)  # Returns Yinterp[1…J] at Xinterp[1…J] which agrees with X[1…N], Y[1…N].
  6. #
  7. #   function Zero(f, xLeft,xRight) # Returns x between Left,Right such that f(x)=0.
  8. #
  9. #   function Minimize(f, x, xLeft,xRight) # Returns minimum value of f(x).
  10. # 
  11. # 
  12. #  This text file explains and gives examples of
  13. #  some of the routines in the file All Library Routines, which
  14. #  should be Opened before trying any of these examples.
  15. #
  16. #################
  17.  
  18. ########
  19. # Here are the xCOD's that are used.
  20.  
  21. xINTERPOLATE
  22.   # xCOD xINTERPOLATE(N:num; X,Y,dy2,w2:RealArray; Ni:num; Xi,Y, Err:num), interpolates y[1…N] of x[1…n] to the points Xi[1…ni]. dy2[1…N],w2[1…N] are work arrays. Output is Yi[1…Ni]. Err: 0=> no error, -1=> two x[j] are identical; an external program.
  23.  
  24. xZERO
  25.   # xCOD xZERO( prog F(n:num; xy:array); x,xLeft,xRight, tolerance, Err:num), Finds a zero of y(x), where F(2, [x,y]) calculates y. Input:guess for x,ragne xLeft,xRight, tolerance for x. Output:x. Err: -1=> x not in Left,Right range.; an external program.
  26.  
  27. xMINIMIZE
  28.   # xCOD xMINIMIZE( prog F(n:num; xy:array); x ,xL,xR, tol, err:num), finds min of y(x), where F(2, [x,y]) does y(x). In: guess for x, xL<x<xR, y(x) < y(xL) or y(xR), tolerance. Out: x. Err:-1=> x not >xL,<xR; -2=> f(x) > f(xL or xR); an external program.
  29.  
  30.  
  31. #############
  32. # Here are the interface routines.
  33.  
  34. # Given arrays y,x such that
  35. #   1) x is ordered, and
  36. #   2) y=f(x) for some function f, then
  37. # Interpolate returns the points Yinterp = f(Xinterp) for
  38. # an array (or single number) Xinterp. 
  39. # (Which need not be the same size as x or y.)
  40. #
  41. # The first work array, w1, actually returns the
  42. # 2nd derivatives of y(x), which are used in this cubic
  43. # spline interpolation.
  44. #
  45.   function Interpolate(X,Y, Xinterp)  # Returns Yinterp[1…J] at Xinterp[1…J] which agrees with X[1…N], Y[1…N].
  46. .   var n, w1,w2, Yinterp, err
  47. .  #  Input:
  48. .  #             x[1…n] ordered array
  49. .  #             y[1…n] array, y=f(x)
  50. .  #             Xinterp[1…J] co-ordinates to interpolate f(x) at.
  51. .  #  Output:
  52. .  #             Interpolate = an array =f( xInterp[1…J]  )
  53. .   begin
  54. .     w1=y; w2=y;   # reserve storage space.
  55. .     n = size(y)
  56. .     Yinterp = Xinterp
  57. .     Ni = size(Xinterp)
  58. .     xINTERPOLATE(n,x,y,w1,w2, Ni, Xinterp,Yinterp, err)
  59. .     if err<>0 then Print("ERROR in INTERP, two identical x values")
  60. .     Interpolate = Yinterp
  61. .   end
  62. .   
  63.  
  64. # Given a one-line NCII text function f
  65. #  (for example:  f(x) = 2 exp( sin(x)) - 3.2 )
  66. # Zero returns x s.t. f(x) = 0.
  67. #
  68.   function Zero(f, xLeft,xRight) # Returns x between Left,Right such that f(x)=0.
  69. .   var Tol, x, err
  70. .  #  Input:
  71. .  #            f = NCII user function like f(x) = 3x-1.
  72. .  #            xLeft, xRight = numbers bracketing x where f(x)=0.
  73. .  #  Output:
  74. .  #             value of x where f(x) = 0.
  75. .   begin
  76. .     fZeroFunc__ = f     # put function into one zzFuncForZero__ uses.
  77. .     Tol = 1e-6;     #desired tolerance of answer
  78. .     x = (xLeft+xRight)/2   # first guess for x
  79. .     xZero(zzFuncForZero__, x, xLeft, xRight, Tol, err)
  80. .     if err<>0 then Print(" ERROR IN ZERO = "+err)
  81. .     Zero = x
  82. .   end
  83. .   
  84.  
  85. #   Passing functions to external routines is tricky. 
  86. # The ONLY allowed format for an NCII routine is 
  87. # one that looks like this, where
  88. # n is the size of the array xy.  
  89. #   This one has xy=[x,y], and calculuates y=f(x).
  90.  procedure zzFuncForZero__(n,xy) # Used by function Zero.
  91. .   begin
  92. .     xy(2) = fZeroFunc__(xy(1))
  93. .   end
  94. .   
  95.  
  96. # Given another one-line NCII text function f
  97. # Minimize returns the minimum value of f(x)
  98. # inside xLeft, xRight.  This particular interface
  99. # to xMINIMIZE only works if the function value 
  100. # at the midpoint of xLeft, xRight is lower than 
  101. # the function at the ends.
  102. # 
  103. # x is changed to the location of the min, and
  104. # the height f(x) of the minimu is returned.
  105. #
  106.   function Minimize(f, x, xLeft,xRight) # Returns minimum value of f(x).
  107. .   var Tol, err
  108. .  #  Input:
  109. .  #             f = NCII user function
  110. .  #             xLeft, xRight = numbers bracketing the minimum
  111. .  #  Output:
  112. .  #              x such that f(x) is the min.
  113. .  #              Minimize = f(x) = minimum y value.
  114. .   begin
  115. .     x = (xLeft+xRight)/2  # start at midpt.
  116. .     fMinFunc__ = f
  117. .     Tol=1e-6
  118. .     xMINIMIZE(zzFuncForMin__, x,xLeft,xRight, Tol, err)
  119. .     if err=-1 then
  120. .       Minimize="ERROR IN MIN: x must be between xLeft,xRight"
  121. .     else if err=-2 then
  122. .        Minimize="ERROR IN MIN: f(x) >= f(xLeft) or f(xRight)"
  123. .     else if err=-3 then
  124. .        Minimize="ERROR IN MIN: iteration maximum hit"
  125. .     else
  126. .        Minimize=f(x)
  127. .      end
  128. .      
  129.  
  130. #   Procedure to caclulate the min function for the xCOD.
  131.  procedure zzFuncForMin__(n,xy) # Used by function Minimize.
  132. .   begin
  133. .     xy(2) = fMinFunc__(xy(1))
  134. .   end
  135. .   
  136.  
  137.  
  138. ################
  139. # Examples:
  140.  
  141. # First, interpolation.
  142. # Say we have these values of x,y:
  143. x = 0…2 @ .2;  y = cos(x);  
  144. #  table(x,y)   # Evaluating this without the first "#"  gave the table. 
  145. # x                  y
  146. # - - - - - - - - - - - - - - - 
  147.   0                  1
  148.   0.2               0.98007
  149.   0.4               0.92106
  150.   0.6               0.82534
  151.   0.8               0.69671
  152.   1                  0.5403
  153.   1.2               0.36236
  154.   1.4               0.16997
  155.   1.6               -0.0292
  156.   1.8               -0.2272
  157.   2                  -0.41615
  158.  
  159. # To find the values of this function at other x values,
  160. # (pretending for a moment that we don't know it's a cosine)
  161. xOther = [.15, .9, 1.7];    # say we want y(x) at these values.
  162. yOther = Interpolate(x,y, xOther)
  163.   yOther[1…3] = [0.98877, 0.62161, -0.12884] 
  164.  
  165. # The real values of the function should be 
  166. cos(xOther)
  167.   [0.98877, 0.62161, -0.12884] 
  168. # which are pretty durn close.
  169.  
  170.  
  171. # Now, for the zero of a function.
  172. pi = Zero(sin, 3, 3.2)   # find where sin is zero. (This is really at π).
  173.   pi = 3.14159272 
  174. # This is pretty close to the real value:
  175. (π-pi)
  176.   -6.5177165e-8    # This is the error; Zero is set for 1e-6 tolerance.
  177.  
  178.  
  179.  
  180. ########
  181. # Second: find where sin(x)=cos(x)^2.
  182. f(x) = sin(x) - cos(x)^2;      # Define the function we want to examine.
  183. # First we need to find a range where this goes through zero.
  184. x = 0…π@.2; y = f(x);
  185. # table(x,y)   # Evaluating the "Table(x,y)" gives this table:
  186. # x                  y
  187. # - - - - - - - - - - - - - - - 
  188.   0                  -1
  189.   0.2               -0.76186
  190.   0.4               -0.45894
  191.   0.6               -0.11654
  192.   0.8               0.23196
  193.   1                  0.54954
  194.   1.2               0.80074
  195.   1.4               0.95656
  196.   1.6               0.99872
  197.   1.8               0.92223
  198.   2                  0.73612
  199.   2.2               0.46216
  200.   2.4               0.13171
  201.   2.6               -0.21876
  202.   2.8               -0.55279
  203.   3                  -0.83897
  204. # From which we can see that the function is zero between
  205. # .6 and .8, as well as between 2.4 and 2.6. Of course, we could
  206. # have done this with the graph, too.
  207. #
  208. # Then
  209. alpha = Zero(f, .6, .8)
  210.   alpha = 0.66623943 
  211. # Try it:
  212. sin(alpha)
  213.   0.61803399 
  214. [cos(alpha)]^2
  215.   0.61803399 
  216. # They look pretty close…
  217.  
  218.  
  219.  
  220. ############
  221. # Finally, the minimum of a function.
  222. # Usually these should be graphed to find values xLeft, xRight 
  223. # which surround the minumum.
  224. f(x) = 3 + (x-1)^2;      # here's the function to minimize.
  225. # This has a min at f(1)=3.  Try 0.8, 1.3 as brackets.
  226. Minimize(f, xLow, 0.8, 1.3)        # Find the min. 
  227.   3      # This is the returned minimum value.
  228. xLow  # and here where the minimum occurs.
  229.   xLow = 0.9999994 
  230. |xLow-1|/xLow            # The fractional error in x is
  231.    5.960594e-7            # just under 1e-6, which was the 
  232.                                     # the tolerance set in Minimize.
  233.  
  234.